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We numerically investigate the expansion of clouds of hard-core bosons in the two-dimensional 
square lattice using a matrix-product-state-based method. This nonequilibrium set-up is induced 
by quenching the trapping potential to zero and our work is specihcally motivated by a recent 
experiment with interacting bosons in an optical lattice [Ronzheimer et al., Phys. Rev. Lett. 110, 
205301 (2013)]. As the anisotropy of the amplitudes Jx and Jy for hopping in different spatial 
directions is varied from the one- to the two-dimensional case, we observe a crossover from a fast 
ballistic expansion in the one-dimensional limit Jx S> Jy to much slower dynamics in the isotropic 
two-dimensional limit Jx = Jy We further study the dynamics on multi-leg ladders and long 
cylinders. For these geometries we compare the expansion of a cloud to the melting of a domain 
wall, which helps us to identify several different regimes of the expansion as a function of time. By 
studying the dependence of expansion velocities on both the anisotropy Jy/Jx and the number of 
legs, we observe that the expansion on two-leg ladders, while similar to the two-dimensional case, 
is slower than on wider ladders. We provide a qualitative explanation for this observation based on 
an analysis of the rung spectrum. 

PACS numbers: 67.85.-d, 05.30.Jp, 37.10.Jk 


I. INTRODUCTION 

Ultracold quantum gases are famous for the possibility 
of realizing many-body Hamiltonians such as the Hub¬ 
bard model, the tunability of interaction strength, and, 
effectively, also dimensionality [T]. This provides access 
to genuine one-dimensional (ID) and two-dimensional 
(2D) physics as well as to the crossover physics between 
these limiting cases. Moreover, time-dependent changes 
of various model parameters can be used to explore the 
nonequilibrium dynamics of many-body systems (see 12 - 
3] for recent reviews). Timely topics that are investigated 
in experiments include the relaxation and thermalization 
dynamics in quantum quenches [SHU], the realization 
of metastable states m 116] . and nonequilibrium mass 
transport dlHIS] and spin transport [50|. Due to the 
availability of powerful analytical and numerical methods 
such as bosonization |2T|, exact solutions for integrable 
systems |22j , or the density matrix renormalization group 
method P51 - E5] , a direct comparison between theoretical 
and experimental results is often possible in the case of 
ID systems laiioiiisiiiii. 

Strongly interacting many-body systems in two spatial 
dimensions, however, pose many of the open problems in 
condensed matter theory and many-body physics, con¬ 
cerning both equilibrium and nonequilibrium properties. 
The reason is related to the lack of reliable numerical ap¬ 
proaches. Exact diagonalization, while supremely flex¬ 
ible, is inherently restricted to small system sizes |26j . 
Nevertheless, smart constructions of truncated basis sets 
by selecting only states from subspaces that are rele¬ 
vant for a given time-evolution problem have given access 
to a number of 2D nonequilibrium problems (see, e.g., 
[IZIIIH])- The truncation of equation of motions for op¬ 


erators provides an alternative approach |29j . which has 
also been applied to quantum quench problems in the 2D 
Fermi-Hubbard model |30| . Quantum Monte Carlo meth¬ 
ods can be applied to systems in arbitrary dimensions in¬ 
cluding nonequilibrium problems (see, e.g., [?TH55] ). but 
suffer, for certain systems and parameter ranges, from 
the sign problem [34] . Dynamical mean-field methods 
become accurate in higher dimensions, yet do not neces¬ 
sarily yield quantitatively correct results in 2D |35| . 

Regarding analytical approaches, we mention just a 
few examples, including solutions of the Boltzmann 
equation m, flow equations |36| . expansions in terms 
of the inverse coordination number m, semiclassical 
approaches [551 155] . or time-dependent mean-field ap¬ 
proaches [4nH42] such as the time-dependent Gutzwiller 
ansatz (see, e.g., mm)- All these methods have pro¬ 
vided valuable insights into aspects of the nonequilibrium 
dynamics in two (or three) dimensions, yet often involve 
approximations. Recently, the application of a nonequi¬ 
librium Green’s function approach to the dynamics in 
the sudden expansion in the 2D Fermi-Hubbard model 
has been explored [55] . 

Although there have been very impressive recent ap¬ 
plications [461448] of the density-matrix renormalization 
group (DMRG) method [35] to 2D systems, the method, 
in general, faces a disadvantageous scaling with system 
size in 2D [351 ]5S] . Tensor-network approaches [531 - f5T] 
that were specihcally designed to capture 2D many-body 
wave functions are an exciting development, with promis¬ 
ing results for the t — J model [53] . A relatively little- 
explored area of research is the time evolution of 2D 
many-body systems in quantum quench problems using 
DMRG-type algorithms [55H57] . 

In this work, we present the application of a recent 
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extension of ID matrix-product state (MPS) algo¬ 
rithms [HSHHO] that is specifically tailored to deal with 
long-range interactions. Such long-range interactions 
arise by mapping even a short-range Hamiltonian on a 
2D lattice to a ID chain for the application of DMRG. 

Recent experiments have started to study the nonequi¬ 
librium dynamics of interacting quantum gases in 2D 
lattices or in the lD-to-2D crossover [mi^EI]. Mo¬ 
tivated by Refs. [Ml [19], we study the sudden expansion 
of hard-core bosons which is the release of a trapped gas 
into a homogeneous optical lattice after quenching the 
trapping potential to zero. The results of Ref. m show 
that strongly interacting bosons in 2D exhibit a much 
slower expansion than their ID counterpart. In the lat¬ 
ter case, the integrability of hard core bosons leads to a 
strictly ballistic and (for the specific initial conditions of 
Ref. [19]) fast expansion that is indistinguishable from 
the one of noninteracting fermions and bosons. In the 
2D case, it is believed that diffusive dynamics sets in and 
virtually inhibits the expansion in the high-density re¬ 
gion, leading to a stable high-density core surrounded by 
ballistically expanding wings m, similar to the behav¬ 
ior of interacting fermions in 2D m- The characteristic 
feature of these diffusivelike expansions in contrast to the 
ballistic case is the emergence of a spherically symmetric 
high-density core, while the ballistic expansion unveils 
the topology of the underlying reciprocal lattice. 

In our work, we investigate this problem for both 2D 
clusters that can expand symmetrically in the x and y di¬ 
rections [see Fig.j^a)] and wide cylinders and ladders [see 
Fig-Bb)]. We use the ratio of hopping matrix elements 
Jx and Jy along the x and y directions as a parameter to 
study the lD-to-2D crossover. For the 2D expansion in 
the isotropic case Jx = Jy, we clearly observe the emer¬ 
gence of a spherically symmetric core, while for small 
values of Jy < Jx and on the accessible time scales, the 
expansion is essentially ID-like. We further compute the 
expansion velocities derived from the time dependence of 
the radius as a function of JyjJx- 

Since we are, in general, able to reach both longer times 
and larger particle numbers in the case of ladders than in 
2D, we present an extensive analysis of multileg ladders 
and cylinders [i.e., ladders with periodic boundary con¬ 
ditions in the (narrow) y direction] with Ly = 2, 3,4 legs 
[see the sketch in Fig. Bb)]. From the analysis of the ex¬ 
pansion in ID systems P-bj. we expect that the short-time 
dynamics is identical to the melting of so-called domain- 
wall states [H2PH1] . in which half of the system is empty 
while the other half contains one particle per site in the 
initial state [see the sketch in Fig.^c)]. The domain-wall 
melting has been attracting considerable attention as a 
nonequilibrium problem in ID spin-^ systems (see, e.g., 
[62l - [6^ 1. Our results show that this similarity between 
the expansion of clusters and the domain-wall melting 
carries over to the transient dynamics on Ly-\eg ladder 
systems, irrespective of boundary conditions. 

A considerable portion of the discussion in both theo¬ 
retical and experimental papers has focused on the ques- 
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FIG. 1. (Color online) Illustration of initial states and ge¬ 
ometries: (a) central block for the 2D expansion; (b) central 
block of size B x Ly-, and (c) domain wall on a cylinder with 
Ty = 4 legs. 


tion of whether there are signatures of diffusive dynam¬ 
ics in the sudden expansion in 2D, in the dimensional 
crossover [HIM], or on coupled chains m- The anal¬ 
ysis of the expansion of fermions in the 2D square lat¬ 
tice starting from an initial state with two particles per 
site (i.e., a fermionic band insulator) suggests that diffu¬ 
sive dynamics is responsible for the slow expansion in the 
high-density regions m This is expected to carry over 
to the bosonic case, yet there only two-leg ladders have 
been thoroughly studied. In linear response, hard-core 
bosons on a two-leg ladder realize a textbook diffusive 
conductor at high temperatures [niiTi]. thus suggesting 
that diffusion may also play a role in the sudden expan¬ 
sion [70]. Curiously, the expansion velocities measured 
numerically for hard-core bosons on a two-leg ladder ex¬ 
hibit a dependence on Jy/Jx that resembles the experi¬ 
mental observations for the true 2D case [laizn]. Here 
we are able to provide a more refined picture. Our anal¬ 
ysis unveils that the sudden expansion becomes faster 
by going from two-leg to three- or four-leg ladders. We 
trace this back to the existence of heavy excitations on 
the two-leg ladder that are defined on a rung of the lad¬ 
der and are inherited from the Jx ^ Jy limit, which 
cannot propagate in first-order tunneling processes in 
Jx/Jy Conversely, the three- and four-leg ladders pos¬ 
sess single-particle-like excitations, which we dub propa¬ 
gating modes, that have a sufhciently low mass to become 
propagating. This picture provides an intuitive under¬ 
standing of the emergence of slow mass transport in the 
sudden expansion in the initial stages of the time evolu¬ 
tion, complementary to the discussion of diffusive versus 
ballistic dynamics. The reasoning is similar to the role 
that doublons play for slowing down mass transport in 
the ID Bose-Hubbard model [B Uni [73HZ5], which has 
also been emphasized in the case of the Fermi-Hubbard 
model [ZS [77]. Our results raise the question as to 
whether the expansion in both directions in 2D and the 
one-directional expansion on wide ladders and cylinders 
will result in the same dependence of expansion veloci¬ 
ties on Jy/Jx for large Ly. It appears that the ladders 
and cylinders, at least for small Ly, preserve some de¬ 
gree of one-dimensionality. A possible scenario is that 
the expansion velocities in the x direction will depend 
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nonmonotonically on Ly for a fixed value of Jy/Jx if ever 
they become identical to the behavior on the 2D systems. 
As a caution, we stress that long expansion times may 
be necessary to fully probe the effect of a 2D expansion 
at small Jy ^ Jx since the bare time scale for charge dy¬ 
namics in the y direction is set by 1/Jy, as pointed out 

in [75]. 

Apart from the nonequilibrium mass transport of 
strongly interacting bosons, there are also predictions 
for the emergence of nonequilibrium condensates at finite 
quasimomenta in the sudden expansion in a 2D square 
lattice. These predictions are based on exact diagonal- 
ization for narrow stripes m, as well as on the time- 
dependent Gutzwiller method [151 jUj. The dynami¬ 
cal condensation phenomenon has first been discussed 
for ID systems (where it actually is a quasicondensa¬ 
tion (HOj), where it was firmly established from exact 
numerical results isni isD and analytical solutions |64] 
(see also [7(71 l82U84] l and has recently been observed in 
an experiment m- In the sudden expansion of hard¬ 
core bosons in ID, the dynamical quasicondensation is 
a transient, yet long-lived phenomenon [TOIISOI as ulti¬ 
mately the quasimomentum distribution function of the 
physical particles approaches the one of the underlying 
noninteracting fermions via the dynamical fermionization 
mechanism [5K115H] . 

It is therefore an exciting question whether a true 
nonequilibrium condensate can be generated in 2D. Our 
results cannot fully clarify this point, yet we do observe 
a bunching of particles at certain nonzero momenta in 
the quasimomentum distribution after releasing the par¬ 
ticles whenever propagating modes as discussed above 
are present. For the melting of domain walls, the occu¬ 
pation of most of these modes, at which a nonequilib¬ 
rium condensation is allowed by energy conservation and 
at which a bunching occurs, saturates at long expansion 
times. The notable exception are certain modes on the 
Ly = A cylinder. This behavior, i.e., the saturation is 
markedly different from the ID case of hard-core bosons 
in the domain-wall melting, where the occupation con¬ 
tinuously increases. The reason for this increase is that 
the semi-infinite, initially filled half of the system will 
indefinitely feed the quasicondensates [HIM]. As such 
an increase is a necessary condition for condensation, we 
interpret the saturation of occupations as an indication 
that either breaking the integrability of strictly ID hard¬ 
core bosons or the larger phase space for scattering in 2D 
inhibits the dynamical condensation of expanding clouds. 
However, even in those cases on the ladder, in which we 
do not see a saturation, the increase is slower than the 
true ID case, suggesting that coupling chains, in general, 
disfavors condensation. Yet a decisive analysis of this 
problem will require access to larger particle numbers 
and times in numerical simulations or future experiments. 
Note that multileg ladder systems can be readily realized 
with optical lattices, using either superlattices [57] or the 
more recent approach of using a synthetic lattice dimen¬ 
sion [55H5D] . Using a synthetic lattice dimension [55], it 


is in principle possible to obtain cylinders, i.e., periodic 
boundary conditions along the (narrow) y-direction. 

The plan of this paper is the following. In Sec. [IT] we 
introduce the model and definitions. Section [Til] provides 
a discussion and definitions for various measures of ex¬ 
pansion velocities employed throughout our work, while 
Sec. jlVj provides details on our numerical method. We 
present our results for the 2D case in Sec. [Vj while the 
results for multileg ladders and cylinders are contained 
in Sec. We conclude with a summary presented in 
Sec. m while details on the extraction of velocities and 
on the diagonalization of rung Hilbert spaces are con¬ 
tained in two appendixes. 

II. MODEL AND INITIAL CONDITIONS 

We consider hard-core bosons on a square lattice and 
on multileg ladders. The Hamiltonian reads 

^ = - H + h.c.) 

+ +^-c-)] • (1) 

Here aj ^ denotes the creation operator on site i = 
(ix,iy) and Jx{Jy) are the hopping matrix elements in 
the x{y) direction. We choose the hopping matrix el¬ 
ement Jx in the x direction and the lattice constant a 
as units and set h to unity; the ratio Jy/Jx is dimen¬ 
sionless. Note that the Hamiltonian is equivalent to the 
spin-i XX model. In ID {Jy = 0), the Jordan-Wigner 
transformation maps the bosons to free fermions m- Lx 
and Ly denote the number of sites in the x and y direc¬ 
tion, respectively. 

We consider different geometries, namely (i) a small 
square-shaped cluster of Lx = Ly = 12 sites with open 
boundary conditions in both directions, (ii) ladders with 
Lx = 60, Ly S {2,3,4} with open boundary conditions 
(OBCs) in both the x- and y-direction, and (hi) cylin¬ 
ders with Lx = 60, Ly € {2,3,4} with periodic bound¬ 
ary conditions (PBCs) in the y direction and OBCs in 
the X direction. For two-leg ladders, the only difference 
between the Hamiltonian with OPC and PBC along the 
y direction is thus a factor of two in the tunneling matrix 
element Jy. In praxis, we obtain the behavior with PBCs 
by just taking the OBCs data with Jy Jyl‘2.. 

For all simulations, we start the expansion from a prod¬ 
uct state, 

i^o)=’ (2) 

iGB 

in real space. To model the fully 2D expansion, we choose 

to be a square-shaped block oi B y. B sites centered 
in the cluster; see Fig. [^a). On cylinders and ladders, 
we study two different types of B'. (i) a block oi B x Ly 
bosons, centered in the x direction and filling all the sites 
in the y direction as shown in Fig.j^b), and (ii) a domain 
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wall, where the left half of the lattice is occupied by a 
block of Lx!“2. X Ly bosons while the right half is empty; 
see Fig. [^c). 


III. DEFINITIONS OF EXPANSION 
VELOCITIES 

There are several possible ways of defining the spatial 
extension of an expanding cloud and thus also several 
different velocities. 


A. Position of the fastest wave front 

One can define the cloud size from its maximum exten¬ 
sion, i.e., from the position of the (fastest) wave front. 
The velocity derived from this approach will typically 
simply be the fastest possible group velocity (provided 
the corresponding quasimomentum is occupied in the ini¬ 
tial state). Thus, this velocity will not contain informa¬ 
tion about the slower-moving particles and any emergent 
slow and possibly diffusive dynamics in the core region. 
We do not study the wave front in this work. 


C. Core expansion velocity 


In the related experiments with ultracold atoms m 
[19], the focus was on the core expansion velocity that 
is derived from the time evolution of the half width at 
half maximum rdt). The reason is that in these experi¬ 
ments, an average over many ID or 2D systems is mea¬ 
sured. Moreover, the core expansion velocity is primarily 
sensitive to the dynamics in the high-density core (but 
insensitive to the ballistic tails) and thus yields slightly 
different information. In case of multiple local maxima, 
the two outermost points are taken. Since in our simula¬ 
tions we have smaller particle numbers compared to the 
experiments [nun], we use linear splines to interpolate 
the density profile between the lattice sites in order to get 
values for rdt) to a better accuracy than just a single lat¬ 
tice constant. The core expansion velocity is defined as 
the time derivative 


_ drc{t) 
“ dt 


(5) 


The full time dependence of Tc and the extraction of Vc 
is discussed in Appendix [A} 


B. Radial velocity 

Theoretically, it is natural to define the radius R as 
the square root of the second moment of the particle 
distribution {ni{t)). Suppose we are interested in the 
expansion in x direction: We average the density profile 
over the y direction to calculate the radius 

Rlit) = (3) 

ix 

where i^a is the center of mass in the x direction and N 
is the total number of bosons. An analogous expression 
is used to define Ry. To get rid of an initial constant 
part, we use i?^(t) = Ryit) — R^(t = 0) to define the 
radial velocity 


with y, = x,y. The corresponding velocity has contribu¬ 
tions from all occupied quasimomenta. It will ultimately 
be dominated by the fastest expanding particles, and for 
the sudden expansion, Ry will be linear in time in the 
limit in which the gas has become dilute and effectively 
noninteracting. 

The radial expansion velocity of ID systems was stud¬ 
ied for the Fermi-Hubbard model [S2] , the Bose-Hubbard 
model [m [70], and the Lieb-Liniger model [93]. For 
Bethe-integrable ID systems, it can be related to dis¬ 
tributions of rapidities |Mj. For a recent study of the 
radial velocity in the 2D Fermi-Hubbard model, see |3S| . 


IV. NUMERICAL METHOD 

Although the Hamiltonian Eq. Q itself is short 
ranged, long-range interactions arise by mapping the 2D 
lattice to a ID DMRG chain. The presence of such long- 
range interactions renders most of the existing DMRG- 
based algorithms for the time evolution (251 IHHHHn] inef¬ 
ficient because a direct Trotter decomposition of the ex¬ 
ponential is not possible. In our work, we use a recently 
developed extension of an MPS-based time-dependent 
DMRG algorithm that is particularly suited for such sys¬ 
tems [53] . The method is based on a local version of 
a Runge-Kutta step which can be efficiently represented 
by a matrix-product operator (MPO) [SS]. The actual 
time evolution can then be performed using standard al¬ 
gorithms that apply an MPO to a given MPS |2S|. An 
advantage of the method is that it can be easily imple¬ 
mented into an existing MPS based DMRG code and has 
a constant error per site. 

For our simulations, we choose the DMRG chain to 
wind along the y direction in order to keep the range of 
the interactions as small as possible (namely Ly). Sources 
of errors are the discretization in time and the discarded 
weight per truncation of the MPSs after each time step. 
The time steps are chosen small enough to make the er¬ 
ror resulting from the second-order expansion negligible. 
We furthermore choose the truncation error at each step 
to be smaller than 10“^°, which is sufficient to obtain 
all measured observables accurately. The growth of the 
entanglement entropy following the quench requires in¬ 
creasing the bond dimension y with time. Conversely, 
since we restrict the number of states to y < 2000, we 
are naturally limited to a finite maximum time tm at 
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FIG. 2. (Color online) Density profiles for the 2D expansion 
from a 4 X 4 cluster with (a)-(c) Jy/Jx = 0.2 and (d)-(f) 
Jy/Jx = 1.0 at times tJx = 0.0,1.0,1.5. 


FIG. 3. (Color online) Radial velocity Vr,x/y in the x direction 
(top three solid lines) and the y direction (dashed lines) for the 
2D expansion from B x B clusters. The small green triangles 
show the result of an extrapolation to B = oo using Eq. 


which the truncation error becomes significant. Note that 
the bond dimension x required for the simulations grows 
exponentially with time. Increasing the particle num¬ 
bers and Ly leads to a faster growth of the entanglement 
entropy and thus to a shorter maximal time tm- How¬ 
ever, we stress that we clearly reach longer times and 
larger systems than is accessible with exact diagonaliza- 
tion (i.e., pure state propagation using, e.g., Krylov sub¬ 
space methods). 

V. TWO-DIMENSIONAL EXPANSION 
A. Density profiles 

We first characterize the expansion by analyzing the 
time- and position-resolved density profile = 

{fii i (t)), where hi i = a] , a, , is the number oper- 

\ X j y \ / / ' X: y ^xi^y ^xi^y 

ator. We present exemplary density profiles for three dif¬ 
ferent times and and two anisotropies Jy/Jx € {0.2,1} 
in Fig. § For small Jy/Jx = 0.2 [Figs.[^a)-( c)], there is 
a fast expansion in the x direction and nearly no expan¬ 
sion in the y direction. This is expected since the bare 
timescale for the expansion in the y direction set by 1/Jy 
is here much larger than the one in the x direction |78| . 
On the other hand, for Jy = Jx, we find four “beams” of 
faster expanding particles going out along the diagonals. 
These beams are even more pronounced for initial states 
with smaller clusters of 2 x 2 and 3x3 bosons (not shown 
here). 

The most important qualitative difference between the 
density profiles at Jy/Jx = 0.2 and Jy/Jx = 1 is the 
shape. In the former case, the profiles retain a rectangu¬ 
lar form, reflecting the underlying reciprocal lattice and 
the different bare tunneling times in the x versus the y 
direction. For the isotropic case, the initial square shape 
of the cluster changes into a spherically symmetrical form 


in the high-density region. This observation is consistent 
with the experimental results of [19]. 


B. Radial velocity 

In order to compare the expansion for different val¬ 
ues of Jy/Jx more quantitatively, we extract certain in¬ 
tegrated quantities from the profiles, which contain rele¬ 
vant information. One such quantity is the radial velocity 
Vr.x/y derived from the reduced radius Rx/y [see Eq. (^]. 
Details on how we extract from the time-dependent 
reduced radius R{t) can be found in Appendix | A[ 

The radial velocities Vr,x and Vr,y for the 2D expan¬ 
sion are shown in Fig. Unfortunately, our simu¬ 
lations for the 2D lattice are restricted to both very 
short times and small numbers of bosons with block sizes 
B G { 2, 3,4 }. For instance, for 4x4 bosons we reach only 
times tjji « 1.5 J“^. The short times prevent us from a 
reliable extraction of the core expansion velocity, which 
would allow for a direct comparison to the experiment 
[I7l[l9]. The experimental results [19] suggest that, for 
increasing Jy, the core expansion velocity in the x direc¬ 
tion decreases dramatically (see Fig. [^, which has been 
attributed to the breaking of integrability of ID hard-core 
bosons [HED]. 

Our results for the radial velocity Vr show that for the 
smallest block size B = 2, tuning Jy/Jx from 0 to 1 
changes the velocity Vj-^x only gradually while the veloc¬ 
ity in the y direction scales almost linearly with Jy. A 
previous study of the expansion of two-leg ladders also in¬ 
dicated that the core expansion velocity exhibits a much 
stronger dependence on Jy/Jx than the radial expansion 
velocity m- We suspect that this weak dependence may 
additionally result from the small number of bosons con¬ 
sidered in our simulations: Increasing Jy allows a hopping 
in the y direction, which reduces the density and thus the 
effective interaction. In other words, tuning Jy/Jx from 
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FIG. 4. (Color online) Momentum distribution function 
nkx,ky (dimensionless) for the 2D expanding cloud of 4 x 4 
bosons at time t = 1.5 J^T^. The solid green lines show the 
solutions to Eq. ([^. 


0 to 1 increases the effective surface of the initial block to 
include the upper and lower boundaries. From the sur¬ 
face, there is always a fraction of the bosons that escape 
and which effectively do not experience the hard-core in¬ 
teraction. This effect becomes more relevant for smaller 
boson numbers, where the bosons are almost immedi¬ 
ately dilute, feel no effective interaction, and, thus, ex¬ 
pand (nearly) ballistically in both directions. For larger 
block sizes B = 3,4, the ratio of surface to bulk is smaller 
and, therefore, interaction effects become more relevant. 
Indeed, we find for B = 3,4 that tuning Jy/ Jx from 0 to 
1 leads to a significant reduction of Vr^x, most pronounced 
for B = 4. 

Even though we have access to only three values of B, 
it is noteworthy that for all values of Jy/ Jx, Vr,xivr,y) de¬ 
creases (increases) monotonically with B and thus with 
total particle number. This tendency is compatible with 
the behavior of the experiments m performed with 
much larger boson numbers, which motivates us to per¬ 
form an extrapolation to B = oo despite the small num¬ 
ber of bosons. We assume that the finite-size dependence 
is dominated by the surface effects of the initial boundary, 
which scales with B. Therefore, we extract the velocity 
for B = oo from a fit to the form 

const 

‘^r,x/y{^) ~ '^r,x/y{,^ ~ 4 ” ( 6 ) 

at fixed Jy/Jx- The resulting values, which are indicated 
by the small green symbols in Fig. should only be 
considered as rough estimates. 


for the 2D expansion. For a purely ID expansion [Jy = 
0), dynamical quasicondensation occurs ai kx = 

[M [sni isi]. As discussed in Refs. mEg], energy con¬ 
servation restricts the (quasi)condensation to momenta 
at which the single-particle dispersion relation e[kx,ky) 
vanishes since the initial state has zero energy, resulting 
in the emission of bosons with, on average, zero energy 
per particle. For a 2D system, this leads to 

€{kx, ky) = —2Jx cos{kxa) — 2Jy cos[kya) = 0. (8) 

The solutions of this equation are indicated by the solid 
green lines in Fig.[^ We indeed observe an accumulation 
of particles at momenta compatible with Eq. <§■ For 
Jy/Jx = 0.2 [Fig. ®a)], there is almost the same weight 
at any momentunifej, compatible with Eq. (|^. We sus¬ 
pect that this is a relict of the short time t = l.b J~^ = 
0.3 J~^ reached in the simulations: Up to this time there 
was almost no expansion in the y direction; thus, we 
have roughly such that is 

initially independent of ky. Nevertheless, closer inspec¬ 
tion shows slightly more weight at compatible momenta 
with k,, = than at those with k,, = 0 even for small 
Jy [see Fig. ®a)]. This becomes much more pronounced 
for Jy = J;j^see Fig. |^c)]. In this case, the strongest 
peaks are at [kx,ky) = These 

four points correspond to the maximum group velocities 
v[kx,ky) = [2Jy;asm[kxa),2Jyasin[kya)) and, in real 
space, manifest themselves via the four “beams” in the 
density profile shown in Fig. [^f). 

Our results do not serve to clarify whether there ac¬ 
tually is a dynamical condensation at finite momenta in 
2D or not since our initial clusters have too few particles 
in the bulk compared to their surface. The fast ballistic 
propagation of the particles melting away from the sur¬ 
face will only be suppressed once the majority of particles 
is in the bulk initially. If we attribute the outermost par¬ 
ticles to the surface, this would require us to be able to 
simulate at least 7x7 clusters. We believe that the accu¬ 
mulation at finite momenta seen in the quasimomentum 
distribution function is due to these fast particles melt¬ 
ing away from the boundary during the first tunneling 
time. Moreover, we would need to be able to study the 
particle-number dependence of the height of the max¬ 
ima in the quasimomentum distribution function or the 
decay of single-particle correlations over sufficiently long 
distances [80]. 


C. Momentum distribution function 

Figure]^ shows the momentum distribution function 

^-i{k,c(,i,ca-j,ca)+ky(iya-jya)) 

:3y 


VI. CYLINDERS AND LADDERS 

In contrast to the 2D lattice, the ratio of surface to bulk 
is much lower for cylinders and ladders, as we initialize 
the system uniformly in the y direction. Moreover, if 
we tune Jy from 0 to 1, the additional hopping in the 
y direction does not lower the density (and with it the 
effective interaction), as it is the case for the fully 2D 
expansion. We thus expect a weaker dependence of the 
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JyjJx = 0.1 JyjJx = 0.5 JyjJx = 1.0 
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FIG. 5. (Color online) Integrated density profiles 
(dimensionless) for the expansion from a 6 x 3 
clnster on a cylinder with Ly = 3. The green dashed lines 
show the location of the half maximum on the left and right. 


results on the number of bosons. Additionally, we can 
reach larger times than for the fully 2D expansion since 
the range of hopping terms after mapping to the DMRG 
chain is smaller. While we can reach times up to tm ~ 
6 for Ly = 2, we are restricted to times up to tm ~ 

4 for Ly = i and tm ~'A Jx^ for Ly = 4. 



domain-wall melting 
regime 

sudden expansion 
regime 

,Jy indep. 

Jy dependent 










6 ti t 2 time 


FIG. 6. Illustration of the time regimes for the expansion of 
blocks (see the text in Sec. VIA for details). 


ID hard-core bosons is diffusion, sustained by numerical 
studies HU. However, in the sudden expansion, the whole 
cloud expands and it is conceivable that the expansion 
appears to be ballistic because the cloud becomes dilute 
too fast, resulting in mean-free paths being on the order 
of or larger than the cloud size at any time IZO]. 

On the other hand, for larger Jy the block in the center 
does not split at ^ 2 , but a region with a high density 
(“core”) remains in the center. The high-density core is 
clearly established already at intermediate Jy/Jx = 0.5, 
where it still expands slowly. For larger Jy, the spreading 
of this core is continuously suppressed. 


B. Integrated current 


A. Density profile 

Figure shows some typical results for the column den¬ 
sity for the expansion of a block on a cylinder with Ly = 
3. We identify three different time regimes for the expan¬ 
sion of blocks, schematically depicted in Fig.|^ First, the 
evolution during the first tunneling time oc l/Ja, is in¬ 
dependent of Jy'. Since we initialize our system uniformly 
in y direction, in the initial longitudinal hopping, there 
cannot be any dependence on Jy and a finite amount of 
time is required before correlations in the y direction can 
build up. 

Then, in a transient regime 0 < ^2 (where ^2 > 
the melting of the block from either side is equivalent 
to the domain-wall melting [111 US] (compare the sketch 
in Fig. [^. From the two boundaries, two “light cones” 
emerge, consisting of particles outside and holes inside 
the block. Both particles and holes have a maximum 
speed of = 2 J^a. Consequently, the time t 2 := B/AJx 
is the earliest possible time at which the melting arrives 
at the center, such that the density drops below one on 
all sites. Thus, ^2 marks the point in time at which den¬ 
sity profiles obtained from blocks start to differ quanti¬ 
tatively from those of domain walls, dehning the third 
time regime. In the case of a ballistic expansion realized 
for Jy Jx, the density in the center drops strongly at 
t 2 and we can clearly identify two outgoing “jets” as two 
separating maxima in the density profiles; see Fig. [^a). 
To be clear, the expectation for the nature of mass trans¬ 
port in a nonintegrable model such as coupled systems of 


In order to investigate the different time regimes fur¬ 
ther, we consider the number of bosons AN{t) that 
at a time t have left the block B where they were 
initialized. This is equivalent to the particle current 


Jl = integrated 

over time and along the boundary dB of the block, 


^N{t) = ni 
i(^B 


= f ds [jX(s) - ■ 

^0 


(9) 


Here by and bi denote the right and left indices ix of the 
boundary of the initially centered block B. We compare 
AN for the expansion on a two-leg ladder starting from 
either central blocks or domain walls in Fig. ^a). To 
this end we normalize AN by the boundary length \dB\, 
which is simply 2Ly a for the central blocks and Ly a for 
the domain walls. 

For short times t < 0.5 J~^ (i.e., t < ti, see the above), 
all curves in Fig. are independent of Jy. For the quan¬ 
tity AN, the first deviations between domain walls and 
cylinders do not occur at ^2 but at 2^2 = B/2Jx, which 
is exactly the time the fastest holes need to travel once 
completely through the block: By definition, AN is not 
sensitive to the density inside the initial block. For the 
expansion of central blocks, particle conservation gives 
a strict bound AN/\dB\ < B/2a, in which case all the 
bosons have left the initial block. These bounds (equal to 
1.5and 3a“^ for H = 3 and B — 6, respectively) are 
approached in the long-time limit of the ballistic expan¬ 
sion for small Jy/Jx = 0.2, which for H = 6, however, 
happens beyond the times reached in our simulations. 
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FIG. 7. (Color online) Comparison of AN/\dB\ (a) on 
a two-leg ladder for the expansion from central blocks 
(dashed/dotted lines, B x 2 bosons) versus the melting of a 
domain wall (solid lines). The lower panels compare domain 
walls on ladders (dashed lines) to domain walls on cylinders 
(solid lines) for (b) Ly = 4 and (c) Ly — 3. The curves for 
the Z/y = 4 cylinder in (b) with JyjJx = 1,2 are nearly on 
top of each other for t > 1.5 


For the domain walls, AN is not bounded (as long as 
the melting does not reach the boundary of the system) 
and grows for small Jy/Jx as AN oc t linearly in time, 
which, via Eq. (§, corresponds to a nondecaying current 
On the other hand, AN gets almost constant for 
large Jy/Jx for both the domain walls and the blocks. 
This indicates that the expansion is strongly suppressed 
on the two-leg ladder, with a high-density core remain¬ 
ing in the center. We speculate that the regime in which 
AN increases only very slowly is indicative of diffusive 
dynamics, by similarity with m- 


C. Propagating modes: Limit of large Jy ^ J^ 


In order to qualitatively understand the suppression of 
the expansion for certain geometries and specific values 
of Ly, it is very instructive to consider the limit of large 
Jy ^ Jx- We discuss this limit in more detail in Ap¬ 
pendix while here, we provide only the general idea 
and discuss the results. The Hamiltonian Eq. 0 can be 
split up into the hopping on rungs (we denote sites with 
the same index ix as a “rung” for both ladders and cylin¬ 
ders), denoted by proportional to Jy, and the hopping 
terms in the x direction proportional to Jx, collected in 
. Our analysis is based on a diagonalization of , 


which is a block-diagonal product of terms operating on 
single rungs. We view the eigenstates of single rungs 
as “modes,” which can be delocalized by . Since a 
coherent movement of multiple bosons is a higher-order 
process of and thus generally suppressed for large 
Jy/Jx, we focus on modes with a single particle on a 
rung. We then look for modes which are candidates for 
a propagation at finite kx- Importantly, the kinetic en¬ 
ergy Ex oc Jx cannot compensate for a hnite Ey cx Jy for 
Jy ^ Jx- Since we initialize the system in states with 
zero total energy, energy conservation allows only modes 
with Ly = 0 to contribute to the expansion in first-order 
processes in Jx/Jy in time. In general, one could also 
imagine to create pairs of two separate bosons with ex¬ 
actly opposite Ey, summing up to 0. Yet cannot 
create such pairs (see Appendix [ b] for details). 

For smaller Jy , the scaling argument of the energy con¬ 
servation does not hold and additional modes (beginning 
with those of small energy Ey) can be used for the prop¬ 
agation in the the x direction; ultimately, for Jy Jx 
any mode contributes to the expansion already at short 
times. We note that modes with strictly Ey = 0 are 
either present or absent at any value of Jy/Jx- 

Such propagating single-boson modes with Ly = 0 do 
not exist on a two-leg ladder: There are, apart from the 
empty and filled rung, only two states with large energies 
Ly = ±Jy. We argue that precisely this lack of modes 
with Ly = 0 leads to the suppression of the expansion 
with increasing Jy/Jx- It is manifest in Fig. 0a) by 
the fact that AN gets almost constant. Thus, we can 
view the expansion to be inhibited by the existence of 
heavy objects (particles of a large effective mass) that 
can propagate only via higher-order processes. This is 
similar to the reduction of expansion velocities due to 
doublons in the strongly interacting regime of the ID 
Bose-Hubbard model ^ mi [701 [21175]. Another effect 
with very similar physics is self-trapping (see, e.g., [44l 

nanTj). 

Whether propagating modes with Ly = 0 exist or not 
depends not only on Ly but also on the boundary condi¬ 
tions in the y direction. This can serve as a test for our 
reasoning. For Ly = 4, we hnd modes with Ly = 0 on a 
cylinder but not on a ladder (see Appendix [b|. We com¬ 
pare AN for these two geometries directly in Figs, [^b) 
and ^c). For small Jy/Jx = 0.2, the additional coupling 
of the cylinders compared to the ladders has (at least 
on the time scales accessible to us) nearly no influence. 
Yet, for large Jy/Jx, we find not only a quantitative but 
even a qualitative difference: For the Ly = A cylinders, 
AN increases linearly in time, irrespective of how large 
Jy/Jx is. Moreover, the slope is (at t > 1.5 J^T^) roughly 
the same for all Jy/Jx /t 0.5 and does almost not de¬ 
crease with time. Using Eq. ([^, we can relate this to 
the presence of a non-decaying current, which we explain 
in terms of an enhanced occupation at momenta com¬ 
patible with Ly = 0. In contrast, on the four-leg ladder 
there are no propagating modes with Ly = 0; thus, we 
expect no linear increase of AN. Indeed, we find that 
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FIG. 8. (Color online) (a),(c) Core expansion velocities Vc,x and (b),(d) radial velocities Wr,a; versus Jy for the expansion of 
a. 6 X Ly block. The left panels (a) and (b) are obtained on Ly = 2,3,4 cylinders; the right panels (c),(d) are obtained on 
Ly = 2, 3, 4 ladders. The green triangles taken from Ref. [19] show the results of the experiments for the fully 2D expansion 
corresponding to the setup of Sec. 


the currents—i.e., the slopes of IS.N in Fig. [^c)— on the 
four-leg ladder decay in time. Yet, the decay is not as 
extreme as for the two-leg ladder, which we explain by 
the existence of modes with lower energies Ey > 0 than 
on the two-leg ladder. For Ly = 3, it is exactly the other 
way around: There are modes with Ey = 0 on the ladder 
but not on the cylinder. In agreement with this, Fig.j^c) 
shows that the expansion on a three-leg ladder is faster 
than on an Ly = 3 cylinder for large Jy/ Jx = 2. 


D. Expansion velocities 


Figure shows the radial and core velocities for the 
expansion of blocks on cylinders and ladders. We note 
that, while Vc and Vr are nearly independent of Jy/Jx 
in the range Jy/Jx = 0.6,..., 1 for the Ly = 4 cylinder 
[Figs. I^a) and [^b)], the values rdt) and R{t) them¬ 
selves actually do decrease when Jy/Jx is tuned from 0.6 
to 1 (see Figs, [l^ and 13 in the Appendixes), due to dif¬ 
ferent short-time dynamics. Further, for the accessible 
times {tm = 3 for Ly = 4), the density profile out¬ 
side the original block is still completely equivalent to 
the domain-wall melting. Nevertheless, rc{t), by defini¬ 
tion, is also sensitive to the maximum value in the center 
of the block, and R{t) is sensitive to the densities at all 
positions. Thus, the velocities shown in Fig. contain 
valuable and complementary information. 

The two-leg ladder (for which the expansion velocity 
has been studied in Ref. IZSI) shows a behavior similar 
to the experimental data for 2D expansions |19| , namely 


that the core velocity Vc drops down to zero with increas¬ 
ing Jy/Jx- However, by comparing different Ly, we find 
a trend towards a faster expansion when Ly is increased 
at fixed Jy/Jx- This trend is in contrast to the naive 
expectation that wider cylinders should mimic the ID- 
to-2D crossover better. In other words, it demonstrates 
that the two-leg ladder does not capture all the relevant 
physics of the expansion in all directions in the lD-to-2D 
crossover, although it shows the same qualitative depen¬ 
dence of velocities on Jy/Jx as the 2D system studied 
experimentally m- However, we understand this from 
our considerations of the limit Jy ^ Jx in Sec. VIC On 
the Ly = 4 cylinder and the Ly = 3 ladder, there exist 
Ey = 0 modes, and thus a preferred occupation of these 
propagating modes with nonzero ky is possible. More¬ 
over, in those other cases in which there are no modes 
with strictly Ey = 0, there are at least modes with lower 

\Ey\ < Jy. 


E. Momentum distribution function 


The momentum distribution on cylinders start¬ 

ing from 6 X Lj, blocks and at fixed time t = 2.0 J~^ is 
shown in Fig. 1^ At small Jy/Jx = 0.2, we observe a 
bunching of particles at the kx = modes indepen¬ 
dent of ky, similar to the fully 2D expansion at the same 
value of Jy/Jx shown in Fig. 

For Jy = Jx and on the Ly = 3 cylinder, the en¬ 
ergy Ey{ky = ±|y) = Jy can be compensated by ki¬ 
netic energy Ex = —2Jx cos{kxa) in the x direction; com- 
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FIG. 9. (Color online) Momentum distribution function 
rik^^ky (dimensionless) for cylinders with (a),(d) Ly — 2, 
(b),(e) Ly = Z and (c),(f) Ly = 4, starting from a, 6 x Ly 
cluster. Data are shown for time t = 2.0 and (a)-(c) 
Jy/Jx = 0.2 and (d)-(f) Jy/Jx = 1.0 (Note that we have a 
symmetry n-kx,ky = nk^,ky = nk^,-ky)- The black dashed 
lines indicate the flat initial distribution at t = 0. 


pare Eq. Indeed, we find a bunching of particles at 
those momenta in Fig. [^e). The Ey{ky = 0) = —2Jy 
and Ey(ky = ^) = 2Jy mode would yield ^ and 

ky = 0, yet we find a slightly higher weight at smaller 
kx in Fig. |^e). However, we note that all these peaks 
for Jy = Jx in Figs. id) and [^e) are not as high as 
their counterparts for Jy/Jx = 0.2. As we have dis¬ 
cussed above and in Appendix there are no modes 
with Ey = 0 for Ly = 2, 3 on cylinders; hence, the max¬ 
ima in njc^^ky are generally suppressed as we go from small 
to large Jy/Jx for Ly = 2, 3. 

On the Ly = 4 cylinder, we find a bunching of particles 
at {kx,ky) = (^, ^) with roughly the same weight for 
all Jy] compare Figs. [^c) and BlL This is in agreement 
with our considerations of Sec. |VIC| since the modes 
with ky = -^ have Ey = 0. The ky = 0, ^ modes are 
suppressed, similar to the case of Ly = 2, 3. 

The question of whether the bunching of particles 
at certain quasimomenta (that requires the existence of 
propagating modes with energies compatible with those 
quasimomenta) will lead to a true dynamical quasicon¬ 
densation at finite momenta can best be addressed using 
the domain walls as initial stats. Here, we are guided 
by the behavior of ID hard-core bosons: In the sudden 
expansion unmni, the dynamical quasicondensation is a 
transient phenomenon, hence the occupation at A: = ± ^ 
first increases and then slowly decreases as dynamical 


FIG. 10. (Golor online) Time evolution of the peak heights in 
the momentum distribution function for cylinders with (a),(d) 
Ly = 2, (b),(e) Ly = 3, and (c),(f) Ly = 4, starting from a 
domain wall. Data are shown for (a)-(c) Jy/Jx = 0.2 and 
(d)-(f) Jy/Jx = 1.0. 


fermionization sets in izni usi US]. The crossover be¬ 
tween these two regimes—the formation and the decay of 
quasicondensates—is given hy t 2 cc B (see also the dis¬ 
cussion in |16j). For the domain-wall melting, the quasi¬ 
condensates are continuously fed with particles with iden¬ 
tical properties due to the presence of an inhnite reservoir 
and thus the quasicondensation peaks in Uk never decay 
but keep increasing. 


Figure[^shows the time dependence of the occupation 
at the maximum of nkx,ky for the domain-wall melting on 
Ly = 2,3,4 cylinders for (a)-(c) Jy/Jx = 0.2 and (d)- 
(f) Jy/Jx = 1. For Jy/Jx = 0.2 and the accessible time 
windows of the Ly = 3,4 cylinders, the occupation indeed 
increases monotonically in time. On the Ly = 2 cylinder 
in Fig.j^a), the maximum initially increases similar as 
for Ly = 3,4, yet for times t > 3J“^ it saturates and 
even decreases, which suggests that no condensation sets 
in. Note that the time scale at which the saturation 
happens is quite large, as it is set by J“^. This suggests 
that there is no condensation even for very small Jy > 0 
on the Ly = 2 cylinder. 

The behavior for Jy/Jx = 1 is quite different. In al¬ 
most all cases, the occupation at the maximum quickly 
saturates, which suggests that no condensation sets in. 
This observation is consistent with the absence of fast 
propagating modes on the L^ = 2, 3 cylinders. Among 
the data sets shown in Fig. 10 'd)-[I^f), there is one ex¬ 
ception, namely the peak atjkx , k^ = (^, ± ^) on the 
four-leg cylinder, which monotonically increases without 
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FIG. 11. (Color online) Time evolution of the occupation of 
the largest eigenvalue Ao (dimensionless) of the one particle 
density matrix for cylinders with (a) JyjJx = 0.2 and (b) 
Jy/Jx = 1.0, starting from a domain wall. The dotted green 
lines show the results of an ID chain [Jy = 0) for comparison. 

a trend towards saturation. This case is thus the most 
promising candidate for a condensation at Jy = J^. 

F. Occupation of lowest natural orbital 

To investigate the question of condensation in more de¬ 
tail, we look at the maximum occupation Aq of the natu¬ 
ral orbitals [SB]. The natural orbitals are effective single 
particle states defined as the eigenstates of the one parti¬ 
cle density matrix (S; Oj). The corresponding eigenvalues 
sum up to the number of particles and can be interpreted 
as the occupations of the natural orbitals. A true con¬ 
densate requires that Aq becomes macroscopically large. 

The largest occupation Aq for the domain-wall melting 
of cylinders is shown in Fig. [m In the ID case, indi¬ 
cated by the green dotted line, the occupation grows, for 
large times, as Aq « 1.38-\/t (BD]. For Ly = 2 we find 
two degenerate natural orbitals with occupation Aq. For 
Jy/Jx = 0.2, we find an initial growth for all Ly = 2,3,4, 
but for Ly = 2, the occupation saturates and even de¬ 
creases for large times t > 3J“^, similar as for the peaks 
in the momentum distribution function. In fact, the 
peaks in the momentum distribution are directly related 
to the natural orbitals with the largest occupation: For 
Ly = 2 there are two degenerate natural orbitals with 
maximal occupation with ky = 0 and ky = ^, and their 
Fourier transformation is peaked slightly above (below) 
kx = ^ ioT ky = 0 {ky = ^). Similarly, for Ly = 3 
{Ly = 4) there are two natural orbitals with maximal 
occupation with fcy = ±|^ {ky = ±^) and one (two) 


with slightly lower occupation with ^ = 0 {kj,= 0, ^), 
leading to the peak structure of Figs.j^b) and[^c) (with 
peaks only at > 0 for domain-wall initial states). 

For Jy = Jx, shown in Fig. [TT|(b), Aq saturates and 
even decreases for the cylinders ot width Ly = 2, 3, but 
keeps growing monotonically for Ly = A (at least on the 
time scale accessible to us), in accordance with Figs.j^f) 
and [l^f). For Ly = 4, we find only two (degenerate) 
natural orbitals with ky = with peaks ed, kx = 

Yet the maximal occupation Aq is significantly smaller 
than in the ID case and seems to saturate at larger times. 

It is instructive to compare Aq to the number of parti¬ 
cles in the expanding cloud AN shown in Fig. defining 
a condensate fraction Xq/AN. AN increases linearly in 
time in ID; hence, the condensate fraction goes to zero 
with 1/y/i, consistent with the absence of true long-range 
order. In the case of cylinders, we never observe a satura¬ 
tion of Xq/AN to a constant nonzero value, but it keeps 
decreasing as a function of time. Therefore, a true con¬ 
densation is not supported by the existing data on any 
cylinder. Yet the survival of a quasicondensation on the 
cylinders is consistent with our data. 


VII. SUMMARY 

Motivated by recent experiments with ultracold bosons 
in an optical lattice [Ml |T9| , we simulated the sudden ex¬ 
pansion of up to 4 X 4 hard-core bosons in a 2D lattice. 
In the limit Jx ^ Jy, we find a fast expansion (at least 
on the time scale accessible to us), similar to the ID case. 
When Jy is tuned to the isotropic limit Jx = Jy, some 
fraction of the particles remains as a high-density core 
in the center and a spherically symmetric shape emerges. 
This trend is compatible with the observations made in 
the experiment of Ref. [T9|. Unfortunately, our results 
for the 2D expansion are dominated by surface effects 
due to the small boson numbers as, in fact, in our sim¬ 
ulations we have more particles at the boundary of the 
initial block than in the bulk. This prevents us from an¬ 
alyzing the core expansion velocity m, yet the radial 
velocities Vr,x decrease monotonically with the block size 
B at any fixed Jy/Jx- We observe a bunching in the 
momentum distribution function at quasimomenta com¬ 
patible with energy conservation. This bunching could 
signal a dynamical condensation at finite quasimomenta 
as in the ID case, where this dynamical quasiconden¬ 
sation [80] has recently been observed in an experiment 
m- Although we cannot ultimately clarify the question 
of dynamical condensation in 2D with our small clusters, 
we believe that the bunching of particles at certain fi¬ 
nite momenta in the 2D expansion Jy « Jx stems from 
surface effects. 

In order to investigate the dimensional crossover fur¬ 
ther, we studied the expansion on long cylinders and lad¬ 
ders with up to Ly = A legs. Correlations between the 
particles in different legs, which lead to a Jy dependence, 
built up on a very short timescale of about one tunnel- 
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ing time in the longitudinal x direction. Up to a time ^2 
that is proportional to the linear dimension of the initial 
block, the expansion of blocks, restricted to either the left 
or right half of the system, is identical to the domain-wall 
melting. On two-leg ladders, the density in the central 
region becomes very weakly time dependent and almost 
stationary for Jyj ^ 1, even for the domain walls. This 
is reflected by a vanishing or even slightly negative core 
velocity, similar to the observations made in experiments 
[mils!. By considering the limit Jy'^ Jx-, we argue that 
this suppressed expansion on the two-leg ladder for large 
Jy! Jx stems from the fact that there are no modes with 
Ey = 0 on single rungs. For cylinders and ladders with 
larger Ly S { 3,4 }, we generically find a faster expansion 
with higher velocities than in the Ly = 2 case. Addition¬ 
ally, there is a dependence of expansion velocities on the 
boundary conditions in the y direction. For instance, the 
expansion on Ly = 4 cylinder is faster than on a four-leg 
ladder. In agreement with our considerations of the limit 
Jy ^ Jx, this is accompanied by a bunching at preferred 
momenta ky = and kx = and an increasing 
occupation of natural orbitals. Yet our data does not 
support a true condensation on any cylinder. 


Finally, we state the interesting question whether the 
expansion velocities on cylinders or ladders will ever show 
the same dependence on Jy / Jx as the width Ly increases 
compared to the expansion of a 2D block. The obvious 
difference is that we fill the cylinders and ladders com¬ 
pletely in the y direction. Due to symmetry, the expan¬ 
sion on cylinders is restricted to be along the x direction 
and, as such, closer to the ID case, at least for small 
Ly. There can thus be two scenarios: Either, even for 
Ly —>■ oo, the velocities of the cylinders might well be 
above the experimental results or, as Ly increases be¬ 
yond Ly = 4, the velocities at a fixed Jy/Jx will depend 
nonmonotonically on Ly. 


Further insight into these questions, i.e., the depen¬ 
dence on Ly or the question of dynamical condensation 
at finite momenta in dimensions higher than one, could 
be gained from future experiments with access to measur¬ 
ing the radius. This could be accomplished using single¬ 
site resolution techniques; see [MM] for work in this 
direction. 
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FIG. 12. (Color online) Reduced radius R{t) for cylinders 
with (a) Ly = 2, (b) Ly = 3, and (c) Ly = 4, starting from 
a. 6 X Ly cluster. The thick dashed lines show the linear fits 
used to extract the radial velocities Vr, which are shown in 
Fig-ib). 


Appendix A: Extraction of core and radial velocity 


Both velocities Vr = and Vc = are time 

derivatives of quantities which are not strictly linear in 
time. Thus, both Vr and Vc themselves are time de¬ 
pendent. Figure shows the time dependence of the 


reduced radius R{t), while Fig. 13 shows the core ra¬ 
dius Tcit). In the ideal case we would expect them 
to get constant in the long-time limit. Unfortunately, 
our calculations are limited to finite times tm = 6 
for the two-leg ladder, tm ~ for Ly = 3 cylin¬ 

ders/ladders, tm ~ 3 for Ly = 4 cylinders/ladders, 
and just tm ~ 1-5 for the 2D lattice. 


The reduced radii all start as R(t) = \/2tJxa on very 
short time scales t < 0.5 J~^. This is clear as we are 
initially confined to the hopping in the x direction, inde¬ 
pendently of Jy. For very small Jy ^ Jx, the reduced ra¬ 
dius remains linear in time with the velocity Vr = \/2Jxa 
at all times, as expected for a ballistic expansion from 
an initial state with a flat quasimomentum distribution 
function [mill- A Jy dependence may show up on a 
timescale t oc Jy^. For larger Jy the slope Vr reduces 
at intermediate times (in the time range where we can 
observe it) but increases again for large tJx. The lat¬ 
ter can be understood as follows: The outermost parts 
have the strongest contribution to the sum in Eq. (|^, 
and naturally these outer parts have the highest velocity 
2 Jxa (and also reached a low density such that they are 
dilute and thus do not see each other any more). As¬ 
suming a fraction p of the particles to expand with v 













13 



10 

9 

8 

7 

6 

5 

4 

3 

2 


t [Jx 


2 


1 2 
t [Jx^] 



FIG. 13. (Color online) Core radius rc{t) for cylinders with 
(a) Ly = 2, (b) Ly = 3, and (c) Ly = 4, starting from a.6x Ly 
cluster. The thick dashed lines show the linear fits used to 
extract the core velocities Vc, which are shown in Fig. |^a). 


us from extracting the core velocity for the 2D lattice, 
where they are too strong at the times reached in the 
simulations. Yet it seems reasonable to extract for the 
cylinders and ladders by linear fits rc(t) = Uc -t + const in 
the same way as for Vr- While it works quite well for the 
ballistic expansion at Jy Jx and quite large Jy^Jx, 
rc{t) still exhibits a stronger time dependence for inter¬ 
mediate Jy, e.g., Jy ~ O.SJx on the Ly = 2 cylinder. In 
the latter case, some of the bosons expand initially during 
the domain-wall melting and thus the block and Tc grow, 
yet then the expansion is slowed down and the extension 
of the high-density block measured by Tc becomes weakly 
time dependent. 


Appendix B: Limit of large Jy ';$> Jx 

We split the Hamiltonian Q into two parts ac¬ 
cording to H = ’^here Hi = 

-Jy Y^y Y collects the hopping terms 

within the rung ix and Hf^ collects the hopping 
terms between neighboring rungs. 


and the rest (1 — p) to form an inert time-independent 
block in the center (see also the argument given in [75]), 
a straightforward calculation shows that R{t) « yRivt 
at large times. This is also the reason why R{t) does 
not settle to a constant value on the two-leg ladder even 
for large Jy, although the core in the center barely melts 
and l\N becomes only weakly time dependent: There is 
always a nonzero fraction of particles which go out from 
the center. 

We extract the time-independent expansion velocities 
Vr shown in Figs. [^and|^by a linear fit R{t) = Vr-t+const 
in the time interval 2.0 J“^ ^ t < tm, where tm is the 
maximum time reached in the simulations; see the above. 
For the 2D lattice, we reach only = 1.5 J“^; thus, we 
fit only in the interval 1.0 J“^ <t< 1.5 J“^ in this case. 
In Fig. we show error bars resulting from similar fits 
but using only the first or the second half of the time 
interval. 

In the time regime 0 < t < t 2 , the core radius is con¬ 
stant, although the cloud already expands: From both 
edges, the block melts, but the location of the half¬ 
maximum density does not move due to particle-hole 
symmetry. Just when the first holes arrive in the cen¬ 
ter of the block, the global maximum decreases and Vc, 
the half width at half maximum, begins to increase. It 
then exhibits strong initial oscillations. The latter stem, 
on the one hand, from the discreteness of the particles’ 
coordinates on the lattice, which is only partly cured by 
the linear splines used to extract Tq- On the other hand, 
the melting of domain walls in ID happens in quantized 
“charges,” which lead to well-defined structures in the 
density profile [531 HO^l I103j . Those oscillations prevent 


1. Two-leg ladder 


In the following we give an explicit expression for 
HI on a two-leg ladder in terms of the eigenstates 
of Hi and We denote the four eigenstates of Hi 

on rung ix as 


|0) = |vac), |1+) 

|2) = at. 2 «Li |vac), |1”) 


qI,i +«l,2 

V2 

al,i - 4^,2 

V2 


vac) 


vac) 


(Al) 


where |vac) denotes the vacuum on rung ix- The corre¬ 
sponding eigenenergies Ey of Hf are listed in Tab. 


We then express a, , and 
states, plug them into 




?^X “I” 1 


Al 


in terms of these eigen- 
and obtain: 


-HY^^+JJx = | 0 ; 1 +) ( 1 +; 0 | + | 0 ; 1 ") ( 1 -; 0 | 

+ | 2 ; 1 +)( 1 +; 2 | + | 2 ; 1 -)( 1 -; 2 | 

+ |1+;1+)(0;2|-|1-;1-)(0;2| (A2) 

+ | 1 +; 1 +) ( 2 ; 0 |-| 1 -; 1 -)( 2 ; 0 | 

-I- h.c.. 


Here, |q;;/ 3) = |a) ® \I3) with a,P G { 0,1+, 1“, 2} de¬ 
notes the tensorproduct of the eigenstates on run gs ix 
and ix + \. The terms in the first two lines of Eq. (A2) 
correspond to just an exchange of the eigenstates a (3 
between the neighboring sites. Thus we can identify the 
terms of the first line to drive the propagation of single 
bosons on top of the vacuum. The second line can be 
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Ly 

= 4 cylinder 

N 

ky[/;] 

Ey [Jy] 

0; 4 

0 

0 

1; 3 

0 

-2 


0.5 

0 


-0.5 

0 


1 

2 

2 

0 

-2.828 


0.5 

0 


-0.5 

0 


1 

0 


1 

0 


0 

2.828 


Ly = 2 ladder 

N 

ky [f] 

Ey [Jy] 

state 

0 

0 

0 

|0) 

1 

0 

-1 

|1+) 

1 

1 

|1-) 

2 

0 

0 

|2) 


Ly = 4 ladder 

N 

Ey [Jy] 

0; 4 

0 

1; 3 

-1.618 

-0.618 

0.618 

1.618 

2 

-2.236 

-1 

0 

0 

1 

2.236 


TABLE Al. Eigenenergies of a single rung. For a given parti¬ 
cle number, degenerate levels are listed by their multiplicity. 


leg ladder, we can distinguish between terms which just 
exchange the eigenstates of neighboring rungs and terms 
which mix them. As on the two-leg ladder, we asso¬ 
ciate the exchange terms with the propagation of modes. 
Since contains only single-particle hopping, the ex¬ 
change terms appear only between eigenstates with N 
and N + 1 bosons on neighboring rungs. Thus, to first or¬ 
der in Jx/Jy, a, mode of N bosons can propagate “freely” 
only in a background of ± 1 bosons per rung. By def¬ 
inition, all these exchange terms do not change the total 
energy Ey. 

For the mixing terms, there is no restriction on the ini¬ 
tial particle numbers on the neighboring rungs. However, 
obviously preserves the total number of particle, 
thus there are only mixing terms for |... A^, ...) o 

|... Af ± 1; Af' =F 1 ■ • The initial melting of the edge 
thus happens via a cascade of subsequent mixing pro¬ 
cesses. For example, consider 


seen as the propagation of a particle on top of a one- 
particle background, or alternatively, a single hole in the 
background of filled rungs. 

In contrast to the terms of the first two lines, the terms 
in the third and fourth row of Eq. (A2) mix different 
eigenstates. If we imagine to start from a domain wall 
|...; 2; 2; 0; 0;...), those are the terms which “create” the 
single particle modes |1^) at the border of the domain 
wall. Subsequently, we would imagine these modes to 
propagate away to the left as single-hole modes and to the 
right as single-boson modes. Yet, for the two-leg ladder 
all these mixing terms change the total energy Ey from 0 
to either +2Jy or —2Jy. Thus, the creation is only possi¬ 
ble via higher-order processes, which are suppressed with 
increasing Jyj A term such as jl'*'; 1“) (2; 0| would not 
change the total energy Ay, but such a term is not present 
in Eq. (A2) due to the conservation of total momentum 
ky. it would change from fcy = 0 + 0 to fcy = 0 -b 

To summarize, we argue that the Ly = 2 ladder is 
special as it possesses the two extremal modes ky = 0 
and TT with a large energy Ey = ±Jy for one particle on 
a rung. Note that energy conservation for large Jy S> Jx 
does not suppress the propagation of the modes |1^) in 
the vacuum, but the creation of these modes at the edges 
of the initial blocks or a domain wall. As a consequence 
the current decays very rapidly as evidenced in Fig. [^a). 


2. Larger cylinders and ladders 

We turn now to the cylinder and the ladder with 
Ly = 4. The eigenenergies of on a single rung 
are listed in Tab. |A1| Giving an explicit expression for 
on an Ly =4 cylinder or ladder is not possi¬ 
ble here, since it contains too many terms. Nevertheless, 
we examine its structure. Similar to that for the two- 


|...4;4;0;0...)^|...4;3;1;0...)^|...4;2;2;0...) 
^|...3;3;2;0...)^|...3;3;1;1...) . (A3) 


On the cylinder there are states with Ey = 0 for any 
number of bosons per rung (see Tab. Al). This makes 
it plausible that cascades like (A31 are possible without 
changing Ey on the single rungs. Indeed, we find the 
corresponding terms in the expression for E[f^ (not 
given here). The initial edge of a block or domain wall 
can thus gradually melt into states with one particle per 
rung while preserving the energy Ey. This is confirmed 
by a strong peak in the momentum distribution function 
depicted in Fig. I^f). These additional fcy = modes 
with Ey = 0, which are not present in the two-leg lad¬ 
der, explain thus the trend of a faster expansion seen as 
higher velocities in Fig. Moreover, we stress that this 
process is independent of Jy, provided that other modes 
with Ey 0 are suppressed and our picture is applicable. 
Indeed, we find that the velocities in Fig. and currents 
(slopes) in Fig. [^b) are roughly independent of Jy/Jx, 
even for moderate Jy/Jx ^ 0.6. 


On the other hand, on the four-leg ladder, there are 
no states with Ey = 0 for one or three bosons on a rung. 
It is thus immediately clear that there can be no mixing 
terms which preserve Ey on every rung separately. More¬ 
over, we find that there are also no mixing terms which 
create modes with opposite energy starting from Ey = 0 
on both rungs. As a consequence, the domain wall melt¬ 
ing on the four-leg ladder requires higher-order processes, 
similar to the two-leg ladder. However, the necessary in¬ 
termediate energies Ey = ±0.613 x 2Jy are smaller than 
for the two-leg ladder, such that these higher-order pro¬ 
cesses processes are more likely. This is reflected in Fig.[^ 
by higher velocities for the four-leg ladder compared to 
the two-leg ladder. 
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